# AG setup x1 <- rep(3600,0); x2 <- x1; event<-x1; y<-x1 subject <- x1; evno <- x1; startag <- x1; # Wei Lin Weissfeld setup x1wlw<-rep(0,1500*8); x2wlw<-x1wlw; ywlw<-x1wlw; eventwlw<-x1wlw; subjectwlw<-x1wlw; evnowlw<-x1wlw; stratwlw <- x1wlw; j<-1; k<-1 for (i in c(1:1500)) { x1t <- rbinom(1,1,0.5) x2t <- runif(1,-2,2) yt <- rexp(1,exp(-x1t+x2t)) if (yt<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-yt; event[j]<-1 y2t <- rexp(1,exp(-x1t+x2t)) + y[j] startag[j] <- 0; subject[j]<-i; evno[j]<-1 x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-event[j]; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-1 j<-j+1; k<-k+1 if (y2t<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-y2t; event[j]<-1 y3t <- rexp(1,exp(-x1t+x2t)) + y[j] startag[j]<-y[j-1]; subject[j]<-i; evno[j]<-2 x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-event[j]; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-2 j<-j+1; k<-k+1 if (y3t<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-y3t; event[j]<-1 y4t <- rexp(1,exp(-x1t+x2t)) + y[j] startag[j]<- y[j-1] subject[j]<-i; evno[j]<-3 x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-event[j]; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-3 j<-j+1; k<-k+1 if (y4t<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-y4t; event[j]<-1 y5t <- rexp(1,exp(-x1t+x2t)) + y[j] startag[j]<- y[j-1] subject[j]<-i; evno[j]<-4 x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-event[j]; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-4 j<-j+1; k<-k+1 if (y5t<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-y5t; event[j]<-1 y6t <- rexp(1,exp(-x1t+x2t)) + y[j] startag[j]<- y[j-1] subject[j]<-i; evno[j]<-5 x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-event[j]; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-5 j<-j+1; k<-k+1 if (y6t<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-y6t; event[j]<-1 y7t <- rexp(1,exp(-x1t+x2t)) + y[j] startag[j]<- y[j-1] subject[j]<-i; evno[j]<-6 x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-event[j]; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-6 j<-j+1; k<-k+1 if (y7t<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-y7t; event[j]<-1 y8t <- rexp(1,exp(-x1t+x2t)) + y[j] startag[j]<- y[j-1] subject[j]<-i; evno[j]<-7 x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-event[j]; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-7 j<-j+1; k<-k+1 if (y8t<1) { x1[j]<-x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 subject[j]<-i; evno[j]<-8 startag[j]<-y[j-1] x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-y[j]; eventwlw[k]<-0; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-7 k<-k+1; j<-j+1 } # y8t } # y7t else { x1[j]<-x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 startag[j]<- y[j-1]; subject[j]<-i; evno[j]<-7; for (l in c(7:7)) { x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-1; eventwlw[k]<-0; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-l k<-k+1 } j<-j+1 } } # y6t else { x1[j]<-x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 startag[j]<- y[j-1]; subject[j]<-i; evno[j]<-6; for (l in c(6:7)) { x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-1; eventwlw[k]<-0; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-l k<-k+1 } j<-j+1 } } # y5t else { x1[j]<-x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 startag[j]<- y[j-1]; subject[j]<-i; evno[j]<-5; for (l in c(5:7)) { x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-1; eventwlw[k]<-0; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-l k<-k+1 } j<-j+1 } } # y4t else { x1[j]<-x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 startag[j]<- y[j-1]; subject[j]<-i; evno[j]<-4; for (l in c(4:7)) { x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-1; eventwlw[k]<-0; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j]; stratwlw[k]<-l k<-k+1 } j<-j+1 } } # y3t else { x1[j]<-x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 startag[j]<- y[j-1]; subject[j]<-i; evno[j]<-3; for (l in c(3:7)) { x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-1; eventwlw[k]<-0; subjectwlw[k]<-subject[j]; evnowlw[k]<-evno[j-1]; stratwlw[k]<-l k<-k+1 } j<-j+1 } } # y2t else { x1[j]<-x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 startag[j]<- y[j-1]; subject[j]<-i; evno[j]<-2; for (l in c(2:7)) { x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-1; eventwlw[k]<-0; subjectwlw[k]<-subject[j-1]; evnowlw[k]<-evno[j]; stratwlw[k]<-l k<-k+1 } j<-j+1 } } # y1t else { x1[j] <- x1t; x2[j]<-x2t; y[j]<-1; event[j]<-0 subject[j]<-i; evno[j]<-1; startag[j]=0; for (l in c(1:7)) { x1wlw[k]<-x1[j]; x2wlw[k]<- x2[j]; ywlw[k]<-1; eventwlw[k]<-0; subjectwlw[k]<-subject[j]; evnowlw[k]<-0; stratwlw[k]<-l k<-k+1 } j<-j+1 } } nobs <- j-1 x1 <- x1[1:nobs]; x2 <- x2[1:nobs] y <- y[1:nobs] ; event <- event[1:nobs] evno <- evno[1:nobs]; startag <- startag[1:nobs] subject <- subject[1:nobs] stratum <- evno stratum <- stratum*(evno>=1 & evno<=7) + 1*(evno==0) + 7*(evno==8) agdata <- cbind(x1,x2,y,event,evno,startag,subject,stratum) nobswlw <- k-1 x1wlw <- x1wlw[1:nobswlw]; x2wlw <- x2wlw[1:nobswlw] ywlw <- ywlw[1:nobswlw]; eventwlw <- eventwlw[1:nobswlw] subjectwlw <- subjectwlw[1:nobswlw]; stratwlw <- stratwlw[1:nobswlw] wlwdata<- cbind(x1wlw,x2wlw,ywlw,eventwlw,subjectwlw,stratwlw) # AG setup library(survival) afit <- coxph(Surv(startag,y,event) ~ x1+x2 + cluster(subject)) # PWP setup stratum <- evno stratum <- stratum#(evno>=1 && evno<=7) + 1#(evno==0) + 7#(evno==8) cfit <- coxph(Surv(startag,y,event) ~ x1 + x2 + cluster(subject) + strata(stratum)) #WLW setup wfit <- coxph(Surv(ywlw, eventwlw) ~ x1wlw + x2wlw + cluster(subjectwlw) + strata(stratwlw)) AG Data x1 x2 y event evno startag subject stratum [1,] 1 1.8604109 0.008245877 1 1 0.000000000 1 1 [2,] 1 1.8604109 0.229454090 1 2 0.008245877 1 2 [3,] 1 1.8604109 0.232493702 1 3 0.229454090 1 3 [4,] 1 1.8604109 0.747638025 1 4 0.232493702 1 4 [5,] 1 1.8604109 1.000000000 0 5 0.747638025 1 5 [6,] 0 -0.4011769 1.000000000 0 1 0.000000000 2 1 [7,] 0 1.2928146 0.022433669 1 1 0.000000000 3 1 [8,] 0 1.2928146 0.104228941 1 2 0.022433669 3 2 [9,] 0 1.2928146 0.324501222 1 3 0.104228941 3 3 [10,] 0 1.2928146 0.357434136 1 4 0.324501222 3 4 [11,] 0 1.2928146 0.437690684 1 5 0.357434136 3 5 [12,] 0 1.2928146 0.721489897 1 6 0.437690684 3 6 [13,] 0 1.2928146 1.000000000 0 7 0.721489897 3 7 [14,] 1 1.0983389 0.522449642 1 1 0.000000000 4 1 [15,] 1 1.0983389 1.000000000 0 2 0.522449642 4 2 WLW DATA x1wlw x2wlw ywlw eventwlw subjectwlw stratwlw [1,] 1 1.8604109 0.008245877 1 1 1 [2,] 1 1.8604109 0.229454090 1 1 2 [3,] 1 1.8604109 0.232493702 1 1 3 [4,] 1 1.8604109 0.747638025 1 1 4 [5,] 1 1.8604109 1.000000000 0 1 5 [6,] 1 1.8604109 1.000000000 0 1 6 [7,] 1 1.8604109 1.000000000 0 1 7 [8,] 0 -0.4011769 1.000000000 0 2 1 [9,] 0 -0.4011769 1.000000000 0 2 2 [10,] 0 -0.4011769 1.000000000 0 2 3 [11,] 0 -0.4011769 1.000000000 0 2 4 [12,] 0 -0.4011769 1.000000000 0 2 5 [13,] 0 -0.4011769 1.000000000 0 2 6 [14,] 0 -0.4011769 1.000000000 0 2 7 [15,] 0 1.2928146 0.022433669 1 3 1 [16,] 0 1.2928146 0.104228941 1 3 2 [17,] 0 1.2928146 0.324501222 1 3 3 [18,] 0 1.2928146 0.357434136 1 3 4 [19,] 0 1.2928146 0.437690684 1 3 5 [20,] 0 1.2928146 0.721489897 1 3 6 [21,] 0 1.2928146 1.000000000 0 3 7 [22,] 1 1.0983389 0.522449642 1 4 1 [23,] 1 1.0983389 1.000000000 0 4 2 [24,] 1 1.0983389 1.000000000 0 4 3 [25,] 1 1.0983389 1.000000000 0 4 4 [26,] 1 1.0983389 1.000000000 0 4 5 [27,] 1 1.0983389 1.000000000 0 4 6 [28,] 1 1.0983389 1.000000000 0 4 7 table(x1,evno*event) x1 0 1 2 3 4 5 6 7 0 705 473 295 197 142 99 80 56 1 783 315 125 57 20 10 2 1 summary(afit) Call: coxph(formula = Surv(startag, y, event) ~ x1 + x2 + cluster(subject)) n= 3360 coef exp(coef) se(coef) robust se z p x1 -0.953 0.386 0.0513 0.0491 -19.4 0 x2 0.932 2.539 0.0268 0.0258 36.1 0 exp(coef) exp(-coef) lower .95 upper .95 x1 0.386 2.594 0.35 0.425 x2 2.539 0.394 2.41 2.671 Rsquare= 0.455 (max possible= 1 ) Likelihood ratio test= 2039 on 2 df, p=0 Wald test = 1968 on 2 df, p=0 Score (logrank) test = 1882 on 2 df, p=0, Robust = 353 p=0 summary(wfit) Call: coxph(formula = Surv(ywlw, eventwlw) ~ x1wlw + x2wlw + cluster(subjectwlw) + strata(stratwlw)) n= 10545 coef exp(coef) se(coef) robust se z p x1wlw -1.54 0.214 0.0528 0.0808 -19.1 0 x2wlw 1.43 4.190 0.0314 0.0471 30.4 0 exp(coef) exp(-coef) lower .95 upper .95 x1wlw 0.214 4.680 0.182 0.250 x2wlw 4.190 0.239 3.821 4.596 Rsquare= 0.292 (max possible= 0.92 ) Likelihood ratio test= 3647 on 2 df, p=0 Wald test = 1106 on 2 df, p=0 Score (logrank) test = 3251 on 2 df, p=0, Robust = 344 p=0 summary(cfit) Call: coxph(formula = Surv(startag, y, event) ~ x1 + x2 + cluster(subject) + strata(stratum)) n= 3360 coef exp(coef) se(coef) robust se z p x1 -0.995 0.370 0.0553 0.0538 -18.5 0 x2 0.969 2.634 0.0323 0.0324 29.9 0 exp(coef) exp(-coef) lower .95 upper .95 x1 0.370 2.71 0.333 0.411 x2 2.634 0.38 2.472 2.807 Rsquare= 0.313 (max possible= 0.998 ) Likelihood ratio test= 1263 on 2 df, p=0 Wald test = 1144 on 2 df, p=0 Score (logrank) test = 1159 on 2 df, p=0, Robust = 583 p=0 afit2 <- coxph(Surv(startag,y,event) ~ x1 + cluster(subject)) # PWP setup stratum <- evno stratum <- stratum#(evno>=1 && evno<=7) + 1#(evno==0) + 7#(evno==8) cfit2 <- coxph(Surv(startag,y,event) ~ x1 + cluster(subject) + strata(stratum)) #WLW setup wfit2 <- coxph(Surv(ywlw, eventwlw) ~ x1wlw + cluster(subjectwlw) + strata(stratwlw)) summary(afit2) Call: coxph(formula = Surv(startag, y, event) ~ x1 + cluster(subject)) n= 3360 coef exp(coef) se(coef) robust se z p x1 -1.02 0.361 0.0513 0.0707 -14.4 0 exp(coef) exp(-coef) lower .95 upper .95 x1 0.361 2.77 0.314 0.415 Rsquare= 0.123 (max possible= 1 ) Likelihood ratio test= 441 on 1 df, p=0 Wald test = 208 on 1 df, p=0 Score (logrank) test = 430 on 1 df, p=0, Robust = 161 p=0 summary(wfit2) Call: coxph(formula = Surv(ywlw, eventwlw) ~ x1wlw + cluster(subjectwlw) + strata(stratwlw)) n= 10545 coef exp(coef) se(coef) robust se z p x1wlw -1.22 0.296 0.0514 0.083 -14.7 0 exp(coef) exp(-coef) lower .95 upper .95 x1wlw 0.296 3.38 0.251 0.348 Rsquare= 0.058 (max possible= 0.92 ) Likelihood ratio test= 635 on 1 df, p=0 Wald test = 216 on 1 df, p=0 Score (logrank) test = 633 on 1 df, p=0, Robust = 157 p=0 summary(cfit2) Call: coxph(formula = Surv(startag, y, event) ~ x1 + cluster(subject) + strata(stratum)) n= 3360 coef exp(coef) se(coef) robust se z p x1 -0.699 0.497 0.054 0.0539 -13.0 0 exp(coef) exp(-coef) lower .95 upper .95 x1 0.497 2.01 0.447 0.552 Rsquare= 0.052 (max possible= 0.998 ) Likelihood ratio test= 178 on 1 df, p=0 Wald test = 168 on 1 df, p=0 Score (logrank) test = 173 on 1 df, p=0, Robust = 195 p=0 # stratum specific estimates summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject) + strata(stratum) + x1:strata(stratum))) summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject),subset=(stratum==1))) summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject),subset=(stratum==2))) summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject),subset=(stratum==3))) summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject),subset=(stratum==4))) summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject),subset=(stratum==5))) summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject),subset=(stratum==6))) summary(coxph(Surv(startag,y,event) ~ x1 + cluster(subject),subset=(stratum==7))) summary(coxph(Surv(ywlw,eventwlw) ~ x1wlw + cluster(subjectwlw),subset=(stratwlw==1))) summary(coxph(Surv(ywlw,eventwlw) ~ x1wlw + cluster(subjectwlw),subset=(stratwlw==2))) summary(coxph(Surv(ywlw,eventwlw) ~ x1wlw + cluster(subjectwlw),subset=(stratwlw==3))) summary(coxph(Surv(ywlw,eventwlw) ~ x1wlw + cluster(subjectwlw),subset=(stratwlw==4))) summary(coxph(Surv(ywlw,eventwlw) ~ x1wlw + cluster(subjectwlw),subset=(stratwlw==5))) summary(coxph(Surv(ywlw,eventwlw) ~ x1wlw + cluster(subjectwlw),subset=(stratwlw==6))) summary(coxph(Surv(ywlw,eventwlw) ~ x1wlw + cluster(subjectwlw),subset=(stratwlw==7))) Stratum 1 2 3 4 5 6 7 PWP-AG -0.79 -0.48 -0.52 -1.12 -0.39 -1.64 0.90 WLW -0.79 -1.13 -1.46 -2.15 -2.45 -3.83 -4.09